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We calculated, within the Gross-Pitaevskii formalism, the critical number of atoms for Bose- 
Einstein condensates with two-body attractive interactions in cylindrical traps with different fre- 
quency ratios. In particular, by using the trap geometries considered by the JILA group [ Phys. 
Rev. Lett. 86, 4211 (2001)], we show that the theoretical maximum critical numbers are given 
approximately by N c — 0.55(£o/|«|)- Our results also show that, by exchanging the frequencies uj z 
and u)p, the geometry with uo p < uj z favors the condensation of larger number of particles. We also 
simulate the time evolution of the condensate when changing the ground state from a = to a < 
using a 200ms ramp. A conjecture on higher order nonlinear effects is also added in our analysis 
with an experimental proposal to determine its signal and strength. 
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Bose-Einstein condensates with attractive interactions 
have been realized with 7 Li since 1995 by the Rice group 
[jlj culminating with experiments that have direct obser- 
vation of the grow and collapse of this condensate 0. 
Measurements of the maximum critical number of atoms 
N c in the condensate, in a trap almost spherical, were in 
good agreement with the theoretical predicted numbers, 
within the experimental uncertainties. 

Recently it has been achieved Bose-Einstein conden- 
sation with 85 Rb H by means of Feshbach resonance, 
which allowed wide tunning of the scattering length, a, 
from negative to positive. The ability to control the scat- 
tering length is used to control and measure the stabil- 
ity condition with the corresponding critical number of 
atoms. 

In Ref. Q it was first shown numerically that for at- 
tractive interactions (negative scattering length a) the 
system becomes unstable if a maximum critical number 
of atoms, N c , is achieved. This limit can be stated in a 
convenient expression by 



N c \a\ 



where m is the mass of the particle confined in a trap 
with frequency uo. k is a dimensionless constant, directly 
associated with the critical number of atoms N c . So, by 
using the above assumption of a spherically symmetrical 
trap, several authors [||]6), including us 0, have calcu- 
lated k with a variety of methods. With the precision 
given in Ref. 0, k — 0.5746. In Ref. ||, it was calcu- 
lated the critical number for a nonsymmetrical geometry, 
but in a case that the frequency ratio is not too far from 
the unity (lo z /lo p = 0.72), giving a result for the number 
of atoms almost equal to the spherical one. 

One can also infer from the variational treatment used 
in Ref. H that the constant k depends on the symmetry 



of the trap. Variational estimates were also considered 
m Ref. lf|. So, in cases of nonspherical symmetry, the 
number k will be dependent on the ratios of the trap 
frequencies, with the equation being scaled by some av- 
eraged frequency. As in mostly of the cases considered 
experimentally the spatial symmetry is almost cylindri- 
cal, with the trap frequencies given by oj x « lo v and u> z , 
we assume lo p — lo x — ui y and a geometrical averaged 
frequency given by Co — (w 2 Wp)^ 3 . 



We define 



A = — 



(2) 



such that the trap will have a "pancake-aspect" if A > 1 ; 
and a "cigar-aspect" if A < 1. The spherical symmetry 
is recovered with A = 1. It is convenient to redefine the 
number k given in Eq. (^), showing explicitly its depen- 
dence on A. In this case, the critical number of atoms N c 
is given by 
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Here, in Eq. ([sj), we observe explicitly the dependence 
of N c in relation to A. By exchanging the frequencies u p 
and LJ Z in the trap, we observe that L — > l z 



lr, — > l p and 



A — ► 1/A. The exchange ratio in this case is given by 



fl(A) = JV ' (w * w ' ) 
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i?(A) is the relevant factor that affects the number of par- 
ticles in the condensate, when exchanging the frequencies 
in a cylindrical configuration. In case that k(X) ~ fc(l/A), 
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we can conclude that tu z > u> p results in a larger number 
of particles inside the trap in the critical limit. 

The above considerations and the numerical calcula- 
tions of k(\) that we are communicating are relevant 
to be taken into account in experiments with BEC in 
cylindrical traps with negative a, like the experiments 
that have been performed in JILA with 85 Rb. The 
JILA's group have considered a "cigar-type" symmetry 
in their approach ||Jl^]. They have determined, recently, 
that k = 0.459 ± 0.012 (statistical) ±0.054 (system- 
atic), for a nonspherical trap, where the frequencies were 
17.24 x 17.47 x 6.80 Hz. Using the above notation, we can 
take ujp = ^Lo x ojy — 2ir x 17.35 Hz. So, the correspond- 
ing value of A used in Ref. jllj was u> z /ujp — 6.80/17.35 
= 0.3919. 

Since the JILA trap is nonspherical, it is worthwhile 
to determine numerically the values of k, for different A. 
Our main goal in the present paper is to systematically 
calculate A: (A) in cylindrical symmetry, cither in "pan- 
cake" (A > 1) or "cigar" type (A < 1), in order to verify 
the favorable geometry of the trap to condensate a larger 
number of atoms, when the two-body scattering length 
is negative. As we are going to show, the slight discrep- 
ancy found by the JILA group, when comparing their 
experimental value of k with the theoretical results, can 
partially be explained by the present study. 

For an atomic system with negative scattering length 
and trapped by an external harmonic oscillator (non sym- 
metric, in general), the Bose-Einstein condensate can be 
described by the Gross-Pitaevskii equation: 



-^ 2 + ^(u 2 x x 2 +uy+ujtz 
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The conditions for the validity of this formalism to de- 
scribe atomic systems with negative scattering lengths 
are given in Ref. |Eq] . Deviations due to quantum fluc- 
tuations and tunneling, that occur near the collapsing 
region, were studied in Refs. [6^13]]. As it appears from 
such studies, the decay probability due to quantum tun- 
neling (that will effectively reduce N c ) is negligible, un- 
less N w N c . 

The wave- function, given by 



#(r,i) = exp(-i/i</^)*(f,0), 



(0) 



where p is the chemical potential, is normalized to the 
number of atoms: 



d 3 r\^>{r,t)\ 2 = N. 



(7) 



p 2 = (2mui/h)(x 2 + 



ujp) and con- 



Using cylindrical symmetry (uj x = 
sidering dimcnsionless units [r = tot, 
V 2 )i C 2 — (2md>/fi,)z 2 ], followed by a new scaling of the 
wave- function, 
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we have 
where 



2 r 2 
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Given the Eqs. (0) and (§), we obtain the normaliza- 
tion of $ to a defined reduced number of atoms n: 
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where, in the critical limit, n = n c = 2V2k. Eq.(g) de- 
pends only on the ratio A = (oj z /u) p ): 
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where f3 = p/(Tioj). So, the normalization constant n, 
given by (|l0|), as well as fc, will depend only on A. 

In our calculation of Eq. (^ we employed the relax- 
ation method propagating in the imaginary time and 
renormalizing $ to 2n at every step We searched 

for stable solutions by varying the number n till a critical 
limit n c . No ground-state solutions are possible for n > 
n c . In Fig. 1, we have the corresponding results for the 
chemical potential as a function of N\a\/lo — n/(2v2)- 
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FIG. 1. The chemical potential, p, is given in units of 
Tuui, as a function of N\a\/lo- Results with spherical symme- 
try (A = 1), in dashed line and with x, are compared with 
results using A = 6.80/17.35 (solid line). Dashed line was 
obtained using shooting-Runge-Kutta method, while the x 
and the solid line were obtained by propagation in imaginary 
time. 
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To obtain the results shown in Fig. 1, we first tested 
our code by running the symmetrical case A = 1 (lu z = 
Up) and comparing the results with very precise ones 
that we have previously calculated with the shooting- 
Runge-Kutta algorithm Q. The plot with x marks 
corresponds to the imaginary time propagation method 
while the dashed line plot refers to the shooting- Rungc- 
Kutta method (in both it was used spherical symmetry) . 
One should note that the unstable solutions (back bend- 
ing branch) are not accessible by the time-dependent 
method. The plot with solid line shows our results for 
a cylindrical symmetry, with the JILA parameters given 
in Q, i.e., uj z = 2tt x 6.80Hz uj p = 2vr x 17.35 Hz. In 
this case, k — 0.552 is approximately 4% lower than the 
spherical case. 

For the propagation, we have used the Peaceman- 
Rachford alternating-direction implicit method . The 
time evolution for cylindrical symmetry was performed 
with a code used in |l6| ]. Our discretization was up to 
200 x 200 in p and £ space directions, and up to 50 in 
the variable r (= Qt), with steps of 0.001. We also con- 
sidered pmax and (max ranging from 2 to 10 depending 
on the symmetry. In the extreme nonsymmetric cases (A 
or 1/A >> 1), the results are more sensible to the grid 
spacing and to these maximum values. In these cases, a 
lack of precision can occur in the third decimal digit of 
the results shown in Table I. 

In Table I, we present the numerical results for the crit- 
ical constant A; as a function of the parameter A = u> z /u> p , 
that can be useful to analyze experiments with differ- 
ent cylindrical shapes. Clearly, the optimal value for k 
occurs for spherically symmetric traps (A = 1), as one 
could also infer from the variational calculations given 
in ||. In particular, we determined the values of k 
for the ratios considered in the JILA experiment |§jtj|: 
The theoretical constant, k « 0.55, is about 4% lower 
than the corresponding number with spherical symme- 
try {k s — 0.5746). This can partially explain the small 
disagreement observed in Ref. JnJ when comparing their 
result with theoretical ones. 

By exchanging the frequencies oj p and uj z in a cylin- 
drical symmetry, it is also shown that the "pancake- 
type" symmetry (u) z > uj p ) is preferable (in order to 
obtain a larger N c ) when fc(A) ps fc(l/A). Considering 
the exchange ratio presented in Eq. (|J) and the results 
shown in Table I for k(X), one can verify the optimal 
geometry to increase the critical number of atoms N c 
trapped in a condensate. Analyzing the "pancake-type 
symmetry" , related with the "cigar- type symmetry" con- 
sidered by the JILA group in Ref. |L1|, we note that 
Ai = 17.35/6.80 = 2.5517, and A 2 = 1/Ai = 0.3919. 
As shown in Table I, both cases will give us practically 
the same constant number k 0.55. So, the relevant 
factor that will decide the convenient symmetry to con- 
densate a larger number of atoms is given by Eq. (^) , in 
this case; and this favors the "pancake-type" geometry: 

i?(Ai = 17.35/6.80) « 1.17 . (12) 



The number of atoms in the condensate can be increased 
by a factor of ~ 17%, just by exchanging the geometry 
of the trap. The above factor can be verified experimen- 
tally, as well as other frequency ratios, with the help of 
Table I and the present relations given for k(X) and R(\). 

We should add that other part of the observed discrep- 
ancy in the experimental value of k could be explained 
by an early collapse of the condensate due to a dynam- 
ical chirp in the wave-function when moving the system 
from a > to a < 0. It means that, when changing 
the scattering length from a positive to a negative value, 
the energy minimum with a > is greater than the cor- 
responding energy minimum with a < 0, such that the 
system will collapse at a lower critical number |l7| . 

We simulated the realistic situation with the parame- 
ters given in Ref. |ll| . We depart from the ground state 
with a = and then ramp it to a < in 200ms. In 
Fig. 2(a) we show the time evolution of the mean-square 
radius, p, for different final negative scattering lengths. 
For a final value of N\a\/lo lower or equal to 0.94fc s , the 
system presents collective excitations; for a larger value, 
the system collapses. So, we conclude that the dynamical 
effects can only account for about 2% of the discrepancy 
observed between the experimental and theoretical values 
of k. This result implies that, the total correction due to 
the nonspherical symmetrical trap and due to dynamical 
effects can only account for a diminishing of about 6% in 
the spherical predicted value of k. For comparison, we 
also present in Fig. 2(b) the corresponding instantaneous 
shift from a = to a < 0. 

TABLE I. Numerical solutions for the critical constant k, 
as a function of A = lo z /uj p . k — k s is for spherical symme- 
try. With (*) we indicate the symmetry considered by the 
JILA group; alternatively, with (|), the corresponding "pan- 
cake-type" symmetry. 



A 


k 


k j k s 


0.01 


0.314 


0.547 


0.02 


0.352 


0.613 


0.05 


0.411 


0.716 


0.1 


0.460 


0.801 


0.2 


0.509 


0.886 


0.3 


0.535 


0.931 


(6.80/17.35) (*) 


0.550 


0.957 


0.5 


0.560 


0.975 


2/3 


0.570 


0.992 


1.0 


0.5746 


1.000 


1.5 


0.570 


0.992 


2.0 


0.561 


0.976 


(17.35/6.80)(|) 


0.549 


0.956 


V8 


0.544 


0.946 


3.0 


0.541 


0.941 


4.0 


0.518 


0.902 


5.0 


0.498 


0.867 


10.0 


0.441 


0.767 


20.0 


0.376 


0.655 


50.0 


0.294 


0.511 
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A larger deviation of k is expected in this CcLSG, cLS this 
numerical simulation [shown in Fig. 2(b)] corresponds 
to a larger chirp in the wave function than in the case 
that a is "ramping" slowly in time. We found that, at 
N\a\/lo = 0.9fc s the system have complex higher mode 
nonlinear oscillations; for a larger value of N\a\/lo, it 
collapses. So, even in this case, we can account to a 
maximum of 10% shift in the value of k (including dy- 
namical and nonspherical effects), when comparing with 
the spherical result. 
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FIG. 2. Time evolution of the dimensionless mean-square 
radius (p) of the condensate when changing the ground state 
from a = to a < 0. We have considered a 200ms 
(r = Cut = 16) linear ramp in (a); and an instantaneously shift 
in (b). In (a), the dashed, solid and dot-dashed lines corre- 
spond to the ramping until N\a\/lo = 0.9k s , 0.94fc s , 0.95k s , 
respectively. In (b), the dashed line corresponds to the ramp- 
ing until N\a\/lo = 0.9k 3 ; and the solid line corresponds to 
the ramping until AT|a|/Zo = 0.91fc s . k 3 is the collapse con- 
stant k — iVjaj//o in spherical symmetry. Trap parameters 
were uj p — 2-k x 17.35 Hz and uj z = 2tt x 6.80 Hz. 

As temperature dependence is being ruled out in 
the experimental analysis, another interesting possibility, 
that could explain a larger deviation in the value of the 
constant k, can be attributed to higher order nonlinear 
effects, that in this case are contributing to increase the 
attractive part of the effective nonlinear potential. The 
relevant effect of a real three-body effective interaction, 
given by a quintic term <E>| 4 <I> in the rhs of Eq. (|lT|) , was 
already pointed out in p8fl . If g 3 is positive, there is a pos- 
sibility of two-phases in the condensate |18|. However, in 
case that 53 has the same negative sign as the two-body 
interaction, one can also obtain a relevant contribution 
that may explain a smaller value for the constant k, as 
it is occurring in the present case. In order to obtain 
the missing part of deviation (~ 10 — 15%), we estimated 
numerically that it is enough to have g% ~ —0.03. 

A way to obtain some definitive conclusion about the 
above conjecture of a relevant role of higher order nonlin- 
earity, is open experimentally by examining particularly 



the case a w 0, when the cubic term in the rhs of Eq. fjll]) 
is replaced by a quintic term. A limit in the number of 
particles at this particular value of a is a good indication 
of negative higher order nonlinearity; and, given N c , the 
corresponding strength of the nonlinear interaction (that 
should mainly come from three-body effects) can be es- 
timated. 
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